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Abstract 

In this paper we develop a model of stochastic gene expression, 
which is an extension of the model investigated in the paper [T. 
Lipniacki, P. Paszek, A. Marciniak-Czochra, A.R. Brasier, M. Kim- 
mel, Transcriptional stochasticity in gene expression, J. Theor. Biol. 
238 (2006) 348 — 367]. In our model, stochastic effects still origi¬ 
nate from random fluctuations in gene activity status, but we precede 
mRNA production by the formation of pre-mRNA, which enriches 
classical transcription phase. We obtain a stochastically regulated 
system of ordinary differential equations (ODEs) describing evolution 
of pre-mRNA, mRNA and protein levels. We perform mathematical 
analysis of a long-time behaviour of this stochastic process, identi¬ 
fied as a piece-wise deterministic Markov process (PDMP). We check 
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exact results using numerical simulations for the distributions of all 
three types of particles. Moreover, we investigate the deterministic 
(adiabatic) limit state of the process, when depending on parameters 
it can exhibit two specific types of behavior: bistability and the exis¬ 
tence of the limit cycle. The latter one is not present when only two 
kinds of gene expression products are considered. 

Keywords: Stochastic gene expression, Pre-mRNA, Piece-wise deterministic 
Markov process. Invariant density 


1 Introduction 


Gene expression and its regulation is a very complex process, which takes 
place in the cells of living organisms, especially in eukaryotes 24 . It is 


widely known that this process depends on the behaviour of crucial sub¬ 
stances, called transcription factors (TFs) and chromatin architecture. Our 
investigation is based on the idea of [^, where a simplihed diagram of gene 
expression was presented. It was mentioned there that genes fluctuate ran¬ 
domly between their activity or inactivity status and transcripts are produced 
in bursts. Stochastic effects at the initial stage are very strong compared to 
both the matter production and degradation processes, so we consider the 
noise of Markov-type origin merely at the activation stage. These claims were 
verihed and analysed through the years [4], 11 , 15 , 16 , 27 . The whole 


scheme describes expression of a single gene, assuming it has n copies, but 
further analyse was performed in the case of one copy only. After activation 
of the gene (which is initiated by binding to the promoter region some of 
TFs), mRNA transcription and protein translation phases follow. At hrst, 
mature mRNA is produced in the nucleus, then it is transported from the 
nucleus to the cytoplasm, where the second phase takes place. As a result, 
new proteins are born. 

In the mentioned class of models, not only transcription and translation 
evolution were considered, but also biological degradation of both types of 
the particles: mRNA and protein. All the processes were recognised as con¬ 
tinuous, so the planar system of ordinary linear differential equations were 
used to represent the dynamics of fluctuations in the level of certain type 
particles. Moreover, hrst equation included stochastic “switch” component, 
being responsible for the control of gene activity status. This system has 
been identihed in Im as a Piece-wise Deterministic Markov Process (PDMP), 


2 







introduced by However, after reflection on these results, an important 
question arises: to what extent does the two-stage model hts the current 
state of biological knowledge? Would adding another stage make description 
of the gene expression more precise? Finally, will the problem be much more 
complicated if we add the third stage? In the mentioned work of there 
is a remark that translated mRNA particle must get through some further 
processing before a new, mature protein is formed. Beside that, plenty of 
thematic books; 22 , 35 and publication sources; [^, claiming that at 
least one additional phase, called primary transcript (or pre-mRNA) pro¬ 
cessing should be taken into account. Actually, in eukaryotic genes, after 
the activation signal, the DNA code is transformed into pre-mRNA form of 
transcript. Then, the non-coding sequences (introns) of transcript are cut 
off. This action is combined with other modihcations widely known as RNA 
processing. Only then we get a functional form of mRNA, which is trans¬ 
ferred into the cytoplasm, where during the third phase, translation phase, 
mRNA is decoded into a protein. In short, we consider three-phase model of 
gene expression with three main components, i.e. three variables xi, 0 : 2 , 2:3 
describing evolution of pre-mRNA, mRNA and protein levels. Firstly, we as¬ 
sume that pre-mRNA molecules are produced at the rate Ai 7 (t), where Ai is 
a constant and we introduce a stochastic binary valued function 7 (t) G { 0 , 1 } 
which marks, at time t > 0 , if the gene is in active ( 7 (t) = 1 ) or inactive 
{■yit) = 0) state. This function will be described in detail in Sec. 2.4 The 


mRNA production rate is equal to A 2 Xi{t), where A 2 is a constant and xi(f) 
denotes the number of pre-mRNA molecules at time t. Similarly, the pro¬ 
tein translation takes the place at the rate A 2 Xi{t), where X 2 (t) denotes the 
number of mRNA molecules at time t. Moreover, all three types of particles 
undergo the degradation process. The total lost of pre-mRNA particles is 
given by diXi = d[xi + A 2 X 1 , where the constant d[ is the degradation rate 
of pre-mRNA particles and another constant A 2 is the rate of converting pre- 
mRNA into mRNA particles. It means that di = d[ + A 2 should be treated 
as the total degradation rate of pre-mRNA particles. This concept takes 
into consideration that pre-mRNA is converted to mRNA 21 , in contrast to 


mRNA which serves as a template for mRNA synthesis, but is not degraded 
during the synthesis. Thus, in other cases we use standard description, i.e. 
the constants ^2 and d^ denote, respectively, mRNA and protein degradation 
rates. This expansion of the previous, simplihed diagram of gene expression 
depicted in 20 , is now presented in Fig. We note that the switching be¬ 
tween active and inactive state of the gene depends on the so-called jump 
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rates (activation/inactivation rates). 


Auto-regulation (possible) 



Degeneration Degeneration Degeneration 

or convertion to 

pre-mRNA 


Fig. 1: An extended scheme of (auto-regulated) gene expression. 


Introduction of the third variable to the model means, that its geome¬ 
try moves unavoidably into space. Although, in the last few years some 
PDMP-based biological models were presented, they focus on the applica¬ 
tions in planar systems: 1^, 20 , 25 , 32 . What is worth mentioning, in 


the paper of it was investigated that three-dimensional Lorenz system 
with stochastic switching “admits a robust strange attractor”, but here we 
concentrate on a situation, when jump rates are not necessarily constant and 
also on the convergence in time of the distribution of the process to the equi¬ 
librium distribution. In the proof of the main theorem, we use some results 
concerning asymptotic stability of Markov semigroups. The main idea is that 
we need to check two conditions: irreducilibity of the semigroup and the ex¬ 
istence of some estimation of the semigroup from below which can be checked 
by using Hormander condition for considered process. An alternative (with 
no reference to Markov semigroups) approach concerning long-time behavior 
of PDMP is based on regularity and convergence results from the papers 
and [^, where similar conditions for the convergence of PDMP were devel¬ 
oped. Recently, asked the question, is it possible to describe long-time 
qualitative properties of the process for spatial dynamics? Here we do such 
an analysis with the aim to include the role of primary transcript in the 
processes basic for eukaryotic genes. 

This paper is organised as follows. First, in Section we present idea of 
the model. Then we discuss the deterministic (adiabatic) limit state, which 
can exhibit two specihc types of behavior: bistability and the existence of the 
limit cycle. Later, we present mathematical description of the process, in¬ 
cluding deterministic part and the stochastic component. Having done that, 
we introduce Markov semigroups and we recall how they can be generated 
by PDMP to describe time evolution of the densities of the process. In the 
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first part of Section we formulate the main theorem of this paper, which 
says that the Markov semigroup related to the model is asymptotically sta¬ 
ble. This means that there exists a stationary density and independently on 
the initial distribution the density of the process converges to the stationary 
density as time goes to inhnity. We hnd a set, an “attractor”, on which this 
three-dimensional distribution is concentrated. The second part of Section 
is devoted to stochastic simulations of the process. We show time-dependent 
and mutual dependent behaviour of levels of pre-mRNA, mRNA and pro¬ 
tein. We also approximate the above-mentioned limit stationary density for 
all types of the particles. In the further part we also refer to the deterministic 
state behavior. In Section we sum up the results of our paper and give 
some conclusion remarks. 


2 The model 


2.1 Construction 

Let xi,X 2 -iX^ denote three non-negative variables, which describe time-evol¬ 
ving levels of pre-mRNA, mRNA and protein, respectively. In accordance 
with the current surveys, we consider that the activation and inactivation rate 
functions go(^i) 2 : 2 ,and gi(a:i, 0 : 2 ,can be constant [^, or can de¬ 


pend on the number of the particles of one type, usually the proteins 8 , 20 


Briefly speaking, the gene is activated with the rate qo{xi,X 2 ^x^) and inac¬ 
tivated with the rate qi{xi,X 2 ,x^). The minimal mathematical assumptions 
about go and qi in the case of two variables are discussed by |^. In line with 
the approach of (^, we study evolution of the following system of ODEs 
with a stochastic component: 


''q qo{xi,X2,X3)^ ^ Q ^i(xi,X2,X3) ^ 

dxi . , ^ 

— = Ai 7 (f) - diXi 

= A2X1 - d2X2 
at 

= A3X2 - dsxs, 


( 1 ) 


where A^ and di are positive constants. 

Remark 1. We pay attention to the fact that a standard three-dimensional 


Goodwin model of an oscillatory gene regulation loop 13 was developed in 


5 








a similar manner to ours and it can be interpreted even in the same way 34 . 


However, instead of the presence of the stochastic process 7 (t), Goodwin 
model contains a non-linear term describing the production rate of mRNA. 
The source of nonlinearity is the dependence of this rate from the protein 
level. In our work we take this fact into account by making the intensity 
functions qo and qi possibly dependent from the level of any type of particles, 
especially the proteins. 

If go and gi are constant, we calculate the expected levels of pre-mRNA, 
mRNA and protein in the molecular population: 


E(a:i) 

E(a; 2 ) 

E(x3 ) 


^1% 


di{qQ + gi) ’ 

did 2 {qo gi)’ 

AiA2A3go 

did 2 d‘i{qQ + gi)’ 


despite the fact that these levels oscillate in time (see Sec. 3.2 for details). 
Using standard rescaling techniques known from investigation of the planar 
model in [^, we obtain the system: 


^ qo{xi,X2,X3)^ ^ ^ ^qi(xi,X2,X3) ^ 

dxi 

dX2 , s 

— = 

dxs 

— = b{X2-X3), 


( 2 ) 


a, 6 > 0, in addition a 7 ^ 6 , a 7 ^ 1 and 6 7 ^ 1 . We investigate this system in 
the next sections of the paper. Our results remain true also if a = 6 or a = 1 
or 6 = 1 (see Remark 3.) 

Remark 2. We can consider a more complicated process containing larger 
number of intermediate steps which lead to equations like in the system 
and the approach presented below will not change. However, we analyse 
the system with three equations for the brevity of notation. Larger num¬ 
ber of intermediate steps introduce time delay, which in the case of negative 
feedback makes the system oscillatory. We have observed such oscillatory be¬ 
haviour even in the three-dimensional case but for very special rate functions 


go and gi (see Fig. 12). 
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2.2 The adiabatic limit 


We shall consider particularly interesting behavior of our model, when both 
of the jump rates qo and qi tend to inhnity, unlike their ratio. In this case, we 
can replace the stochastic process 7 (t) by its expected value V ;= Ey = 

to obtain a state called deterministic or adiabatic limit. Hence, the 
system transforms to deterministic system of three ODEs: 



r(a;3) - xi 
a{xi - X2) 
b{x2 - X3). 


(3) 


Depending on the values of the parameters a, b we investigate some specihc 
types of behavior 14 . Firstly, we consider the case of the positive autoregu¬ 
lation, i.e. when E is an increasing function of 0 : 3 . Assume that the equation 
r(c) = c has three roots Ci < C2 < C3 in the interval ( 0 , 1 ) and r'(ci) < 1 , 
r'(c 2 ) > 1, r'(c 3 ) < 1. Then the system (|^ has three stationary points 
Xj = (ci,Ci,Cj), i = 1,2,3. The linearization of ([^ at Xj leads to the follow¬ 
ing characteristic polynomial 


F(A) = (1 + A)(a + A )(6 + A) - T'{ci)ab. 


Since r'(ci) > 1 we have P{0) < 0, and P{X) > 0 for sufficiently large, 
which means that the polynomial P has a positive root and, therefore, the 
point X 2 is unstable. Now we check that stationary points are xi and X 3 are 
asymptotically stable, i.e. all roots of P have negative real parts. Indeed, if 
a = b = 1 then P has all roots with negative real parts: 

Ai = r'(c,)'/'-!, A2 = -l+(-|+4^)^'(c,)'/^ A3 = -i+(-|-^)r'(Q)'/^ 

If we find some coefficients a, b such that P has a root with a nonnegative 
real part, then we find some coefficients a, b such that P has a root with a 
zero real part, i.e., A = ai, a G M. But then a = 0 or 

{a + b + ab) = a^ and (1 -|- a -|- b)a^ = (1 — V'{ci))ab. 

Both cases are impossible if r'(cj) < 1, hence we obtain a bistable state 
(see [^). In such a case one can expect that the stationary density of the 
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stochastic process will be bimodal provided that go and gi are hnite, but suf- 
hciently large. We check this by performing extensive numerical simulations 
presented in Sec. |3.2| In the case of the negative autoregulation, i.e. when 


r is a decreasing creasing function of 0 : 3 , we have only one stationary point 
X = (c, c, c), where c is the unique solution of the equation r(c) = c. Observe 
that in the case a = h = 1 the polynomial P has one negative real root, and 
two complex roots with positive real parts if r'(c) < - 8 . This suggest that 
in this case the limit cycle can appear. We check its existence by simulating 
the system in Sec. 3^ This case is especially interesting, since from the 


standard Bendixson—Dulac theorem it follows that such limit cycle oscilla¬ 
tions are not observed in the two-dimensional system studied before. Again, 
one can expect that for sufficiently large go and gi the stationary density for 
the process will be distributed close to the limit cycle trajectory. 


2.3 Two deterministic systems 

For a hxed state of the gene, which determines the value of 7 (f) = i, i ^ 
{ 0 , 1 }, the process is purely deterministic and we get the system of the hrst 
order differential equations 


" dxi 




= i — Xi 


dt 

dx2 , , 

— = a(x, - xj 

dxs 

n = - "=>)■ 


(4) 


with the initial condition xq = G with a ,6 > 0 , a 7 ^ 6 , a 7 ^ 1 

and 6 7 ^ 1. The solution 7r*(xo) of this system is 


7 r-(xo) = ilexp (Mf)(xo -'?!), (5) 


where 1 = [ 1 , 1 , 1 ] and 


M = 


-10 0 
a —a 0 
0 6—6 


Moreover, with a similarity to the two-dimensional case [^, we have: 


7(xo) = 1 -7(1) + 7(xo). 


( 6 ) 








In Fig. phase portraits of the system Q for both values of i G {0,1} 
are shown. Each time, there exists one stationary solution: for z = 0; a point 
( 0 , 0 , 0 ) is asymptotically stable steady state, as is a point ( 1 , 1 , 1 ) for i = 1. 
Looking at the right-hand sides of system]^ we state that if Xi > 1, then (no 
matter what the value of i is), Xi decreases, moreover X 2 as well x^ follow 
Xi- On the other hand, if Xi < 1, then it stays in the interval [0,1] forever, 
oscillating between 0 and 1, the same happens with X 2 and X 3 . Hence, we 
can reduce the phase space of the process to a cube X = [0,1]^. 



Fig. 2: A sample solutions of Eq. (j^ for a = 2 ,6 = 3, z = 0 (left) and 
a = 2,6 = 3,z = l (right). 


2.4 PDMP: a definition 

We will briefly mention an idea behind PDMP introduced in 0 - We 
consider qQ{xi,X 2 ,x^) and qi{xi,X 2 ,x^) as two continuous and non-negative 
functions on such that: 

go ( 0 , 0 , 0 ) ^ 0 and gi ( 1 , 1 , 1 )^ 0 . 
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Let zq G {0,1}, To = 0, xq G and we define a (random) function 
7 : [ 0 , cxd) — >■ { 0 , 1 } satisfying 7 ( 0 ) = zq and 


7(f) : = 


if Tn-i ^ f < T„, 


1 - z, if f = Tn, 
where for n ^ 1 , T„ is a positive random variable satisfying: 


( 7 ) 


X, := 


where 


: Prob(T„ - Tn-i ^ t\ 7(7k-i) = i) 


f f' ^ 

(8) 

: 1 - exp - / qi{7ii{s,:s.s))ds , 

\ Jo J 


if Tn-1 ^ S <Tn, 

M if S = Tn, 

( 9 ) 

X 

3 

1 

1 

3 

1 

(10) 


In consequence, replacing a constant value z G {0,1} in the system (|^ by a 
stochastic process 7 (f): 


( dxi , , 

dx2 , . 

— = .(x, - xj 

dxs 

^=6(^2-X3). 


( 11 ) 


gives a definition of a Markov process (C(f) called a piece-wise deterministic 
Markov process, described by the quartet: 


C(f) := {xi{t),X 2 {t),X:i{t),^{t)) = (x(f), 7 (f)). (12) 

The state space of this process is X = X x {0,1}. The remaining character¬ 

istics are the jump rates g* and the jump distribution J((x, z), •) being the 
Dirac measure ^(x,i-i) such that 

J((a;, z),X) = 1. (13) 

A random variable is called a time of the n-th jump of the process. In 
it was shown that in such a case = T^. — T^-i > 0 , where k ^ 1, < cx) 

and 

lim Tfc = CX), (14) 

k^oo 

which means that the process is well-defined for all times f > 0 . 
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2.5 Markov semigroups and their link with PDMP 


Now we will recall some definitions abont Markov semigronps. We nse them 
to describe the evolution of distributions of the process given by the system 
(|^. Detailed information about some connections between semigroup theory 
and stochastic processes can be found in 19 or [^. Let (X, S,m) be a 
a—finite measure space and let D C = L^(X, E,m) be the set of the 
densities, i.e. 

D = {feL^:f^0, 11/11 = 1 }. 


Definition 1 A linear D preserving mapping P ■. ^ is called a Markov 

(or stochastic) operator. 

Definition 2 A family {P(t)}t^o of Markov operators, which satisfies the 
following conditions: 

• D(0) =Id (identity condition), 

• P{t + s) = P{t)P{s) for s, t ^ 0 (semigroup condition), 

• for each f G the function t —)■ P(t)f is continuous with respect to 
the norm (strong continuity), 

is called a Markov semigroup. 


Definition 3 A Markov semigroup {P{t))t^Q is partially integral if there 
exist to > 0 and a measurable function /c : X x X —)■ M+, such that for every 
feD: 


k{p, q)m{dp)m{dq) > 0 


(15) 


and 


PMfip) S j 


(16) 


Definition 4 A Markov semigroup {P{t)}t^o is asymptotically stable if 

• there exists an invariant density for {P(t)}t^o, i.e. f* E D such that 
Pit)r = r for all t>0, 

• for every density f E D : 
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0 . 


(17) 


lim \ \P{t)f 

t^OO 


r\\ = 


Below we define a property, which is in some sense “opposite” to asymptotic 


stability, introduced in 18 


Definition 5 A Markov semigroup is sweeping (or zero-type) with respect to 
a set A E if for every f E D : 


lim / P{t)f{x)m{dx) = 0. (18) 


A precise instruction on how to construct Markov semigroup for PDMP 
is given by [^. Using the analogy with the two-dimensional model, we write 
Fokker-Planck system of equations for the partial densities /o, fi of the pro¬ 
cess 

I ^ “ ^2)/o) + {{X2 - X3)/o) = qifi - qofo 

(19) 

where /o, fi are the functions dehned on [0, cx)) x [0,1]^ such that for any 
Borel set C M+ x M+ x M+ 


Prob(a;(f) G fB, 7 (f) = i) = 



fi{t, Xi, X 2 , X 3 ) dxi dx 2 dxs, z = 0 , 1 . ( 20 ) 




For the reason of the presence of three spatial variables and a wide range of 
possible jump rates, system (19) is difficult to be solved analytically. How¬ 
ever, we will use Markov semigroup {P(f)t>o} generated by this process to 
prove that it has stationary density, which is an equilibrium with respect to 
time evolution of the distributions. 


3 Results 

3.1 Asymptotic stability 

In this section we present main result of this paper. We consider two par¬ 
ticular solutions of the system (|^. The hrst, 0(f) = (0i(f),02(f),03(i)) is 
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the solution of (|^ with i = 0 and the initial condition (0i(O), 02 ( 0 ), 03(0)) = 
(1,1,1). The second 0(t) = (0i(t),02(^),'03(^)) is the solution of (|^ with 
i = 1 and the initial condition (0i(O),02(0),03(0)) = (0,0,0). We conclude 
that 0 and 0 are two solutions of the system ([^ with i = 0 and i = 1 , 
respectively, which join the asymptotically stable points ( 0 , 0 , 0 ) and ( 1 , 1 , 1 ) 
We construct the set A in the following way. Let Aq be the surface made of 
all solutions of the system ([^ with i = 1, which start from any point lying on 
0. This is also the case with Ai, being the surface made of all the solutions 
of the system (|^ with i = 0, which start from any point lying on 0. We 



Fig. 3: The boundaries of A, Aq - hlled and Ai - transparent for a = 2 and 
b = 10 . 

derive algebraic formulas describing Aq as well as Ai in Appendix A. Having 
done that, we dehne A as a subset of [0,1]^, bounded by Aq and Ai. In Fig. 
we show geometric visualisation of A. For comparison, in Fig. we portray 
the sketch of this set, obtained by numerical simulations of the trajectories 
of the process (see Sec. 3.2). 

Now we can formulate the main result of the paper. 
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Theorem 1 Let ^^(x) >0, i = 0,1 for x G [0,1]^. Then, the Markov semi¬ 
group {Pit)}t^o is asymptotically stable and the support of the invariant den¬ 


sity is the set A = A X {0,1}, where A is expressed in the basis (22) by 
A = {{x-y + z, - y^ + z^) : I ^ x ^ y ^ z ^ Q}. (21) 


A general idea beyond the strict proof of this theorem is provided by . 
However, the proof in the two-dimensional model is simpler, because all the 
properties of the attractor are easy to deduce using geometrical arguments. 
For example, in this case the proof that A is an invariant set for the process 
follows immediately from the Muller theorem and the communication 
between states inside A follows from the Darboux property. Since the geo¬ 
metric arguments in the three-dimensional model are not obvious, we need to 
use precise formula to dehne the set A and to prove its properties. Namely, 
we follow 


31 and prove that A is a set such that 


• A is invariant for the process, i.e. if (a;i(0), 0 : 2 ( 0 ), 0 : 3 ( 0 ), 7 ( 0 )) G A, then 
(xi(f),X 2 (t),X 3 (f), 7 (t)) G A for any f > 0, 

• trajectories (xi(t),X 2 (t),X 3 (f), 7 (f)) of the process starting from any 
arbitrary point from [0,1]^ x {0,1} converge to A when time goes to 
inhnity, 

• there is no smaller set satisfying these two conditions above. 

Detailed mathematical proofs of these claims are long and provided in Ap¬ 
pendix B. In Fig. I^we show two-dimensional projections of A onto the 2D 
plane, looking exactly the same as the set proposed by (^. 


3.2 Stochastic simulations 


Although it is difficult to solve Fokker-Planck equations (19) analytically. 


here we discuss stochastic simulations of the trajectories and distributions 
of the system 0 . made to check the accuracy of our statements. Such 


an approach, based on the 12 algorithm, was used by 25 to visualize the 


evolution of the trajectories in the model of self-renewal cells differentiation. 
For our model a similar code in Wolfram Mathematica environment was 
generated and run. 


We have compared the trajectories of the system (11) for selected values of 


the parameters a, b and jump rates qo{xi, X2, X3), gfyxi, X2, X3) up to the hnal 
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cre-mRNA pre-mRNA 

Fig. 4: Projections of A onto the plane, a = 2,b = 3. Left: pre-mRNA and 
mRNA. Right: pre-mRNA and protein. 


time moment T = 150 or at least 1300 jumps were performed. Afterwards, 
we depicted time evolution of the levels of all three kinds of gene expression 
products for a = ^,b = ^,qo = 3 and qi = 6 with the initial condition 
x(0) = ?/(0) = ;2(0) =^, separately in Fig. [^and pairwise in Fig. (left). 
In Fig. I^and in Fig. (right), where we use the same parameters as in the 
previous case, we assume that the inactivation rate depends on the protein 
level, i.e. qi{xi,X 2 ,X 3 ) = qix^. 

We notice that in both cases, as it was expected, fluctuation in pre- 
mRNA level is much stronger than it is in mRNA or protein level. However, 
when the jump rates are constant, then all three levels seem to vary in a more 
limited range than in the case when the jump rates are protein-mediated: the 
values of standard deviation for consecutive phases are 0.151,0.0887,0.0609 
and 0.182,0.099,0.0617, respectively. Moreover, we empirically calculated 
correlation level between each two phases. While pre-mRNA and mRNA 
levels (with the values of the coefficient 0.555 for constant and 0.573 for 
protein-mediated jump rates), as well as mRNA and protein (0.685 and 0.668) 
levels were significantly correlated, pre-mRNA and protein levels were poorly 
related to each other (0.07 and 0.058). 
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Fig. 5: Trajectories of the stochastic process ( |ll[ ). The initial condition is 
a:(0) = 2/(0) = ^,z{0) = 7(0) = 0 and a = = |, 5 'o = 3, gi = 6 are set 

to show the level of pre-mRNA (top), mRNA (center) and protein (bottom). 
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Fig. 6 : Trajectories of the stochastic process ( |ll[ ). The initial condition is 
x(0) = = \,z{Q) = 1 , 7 ( 0 ) = 0 and a = 2 , & = |, go = 3, gi = Q z(t) 

are set to show the level of pre-mRNA (top), mRNA (center) and protein 
(bottom). 
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Fig. 7: Two-dimensional trajectories of the stochastic process (11). The 
initial condition is x( 0 ) = |,i/( 0 ) = ^,z{0) = |, 7 ( 0 ) = 0 and a = ^,b = 
, go = 3 are set to show the dependence between pairs of the variables, for 
constant (gi = 6 , left) and protein-mediated (gi = 6 x 3 , right) jump rates. 
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Fig. 8 : The sketch of A obtained by numerically portrayed trajectories of 
the process in the three-dimensional space. The initial condition is a;(0) = 
|,l/(0) = ^,.2(0) = 1,7(0) = 0 and a = 2 , 6 = 10 for constant (left) and 
protein-mediated (right) jump rates. 


Leaving all parameters unchanged, we analysed the distributions obtained 
by the simulations of system ( 11 ) with constant and linearly dependent 
inactivation rate function gi(xi, X 2 , X 3 ), respectively; see Fig. To fol¬ 
low the behaviour of a gene, we pictured two-phase marginal distributions, 
i.e. p(t, xi,X2), p(t, X2,X3), p(t, xi,X3), where p(t,Xk,Xj) = fo(t,Xk,Xj) + 
fi{t,Xk,Xj). The graphs were made by simulating the system (11) up to 
Tp = 15, repeated 5000 times and hence it is expected (see also | 20 |) that the 
points {x(Tp), y(Tp), z(Tp),i) obtained this way approximate stationary dis¬ 
tributions fi{xi,X 2 ,X 3 ). Fluctuation strength, described by the jump rates, 
decides about the broadness of p: stronger fluctuations give the broader dis¬ 
tribution. However, on the left hand side of Fig. |^the jump inactivation rate 
is twice as large as the jump activation rate, so the inactive state becomes 
dominant. As a result, the distribution much more signihcantly points into 
the zero direction. 

In Sec. 2^ we discussed and justified two specific types of behavior of the 
process in the deterministic (adiabatic) limit. In the first case it is expected 
that when go and gi are sufficiently large and E 7 is an increasing function of 
the protein level X 3 , we should obtain that the stationary density is bimodal, 
see Fig. 10 On the other hand, when both of the jump rates are still 


large, but this time E 7 is a decreasing function of the protein level X 3 , we 
observe the existence of the limit cycle in the deterministic limit. However, 
after simulating the process for a sufficiently long time, we should obtain its 
density distributed close to this limit cycle. A comparison of the distribution 
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Fig. 9: Marginal distributions p{t,Xk,Xj) = fo{t,Xk,Xj) + fi{t,Xk,Xj) calcu¬ 
lated for time t = 15 with simulations of the system ( p!T| repeated 5000 times 
for constant (left) and protein-mediated (right) jump rates. 


0.4 

pte mRNA ® ® 


12.5 

7.5 
5.0 

2.6 
0 


mRNA 


pte mRNA 


20 








Fig. 10: Marginal bimodal distributions p(t,Xk,Xj) = fo(t,Xk,Xj) + 
fi(t,Xk,Xj) calculated for t = 15 with simulations of the system (11) with 
a = b = 1 and jump rate functions qo{x^) = 10(0.01 + gi = 10 ■ 0, 2. 
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Fig. 11: Marginal distributions p(t,Xk,Xj) = fo(t,Xk,Xj) + fi(t,Xk,Xj) cal¬ 
culated for t = 15 with simulations of the system (0 with a = b = 1 and 
jump rate functionsgo = n, gii^Xs) = 9 ■ lO^^nXg^ with n = 10. 
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Fig. 12: The limit cycle for a = b = 1 and jump rate functions Qq = n and 
= 9 ■ 10^*^ n Xg® in the adiabatic limit. 


of the stochastic process (Fig. 0 and its deterministic approximation (Fig. 
12) is provided. The latter case pays particular attention, because it is a kind 


of behavior which is not present in the two-dimensional system from 20 


4 Conclusion 


We have studied a model of stochastic gene expression with the contribu¬ 
tion of three main phases. Our investigation is based on the two-dimensional 


model introduced by 20 , including: activation of the gene, mRNA transcrip¬ 


tion and protein translation. Activity of the gene is regulated stochastically, 
namely by a Piece-wise Deterministic Markov Process, [^. In eukaryotes, 
where the transcript is produced in bursts, we can neglect other sources of 
stochasticity. However, many reports: 1^, 22 , 1^, 36 suggest that at least 


one additional phase, i.e. pre-mRNA level regulation should be considered 
as well. This moves the state space of deterministic part of the process into 
We have analysed long-time behaviour of densities of the process. Using 
Markov semigroup techniques, we have shown that its distribution converges 
to equilibrium, i.e. there exists a stationary density, such that independently 
on the initial distribution its evolution is being stabilised with this density, 
when time goes to infinity. Moreover, we have found a set, an “attractor”, 
which is a support for the equilibrium. Using statistical approach, we visu¬ 
alised the trajectories of the process and approximated stationary distribu- 
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tions in the cases of constant and protein-mediated jump rates. Moreover, we 
discussed two specific types of behavior of the process: bistability and the 
existence of limit cycle trajectory. To summarize, we obtained qualitative 
and statistical results for the long-time evolution of three-phase kinetics of 
the eukaryotic gene. We notice that the main result is in agreement with the 
one from the two-dimensional model. This suggests that a sequence of gene 
transformations described by equations including stochastic activation and 
only production and degradation processes, does not have a significant influ¬ 
ence on stabilizing long-time distribution of the product levels. However, our 
analysis has also shown that it is not entirely true that the three-dimensional 
system has necessarily analogous dynamics to the two-dimensional one stud¬ 
ied before: the appearance of the limit cycle in our system, is not possible in 
the planar model. Thus, larger number of intermediate steps in the case of 
negative feedback can make the system oscillatory. 

Nonetheless, another question is how to formulate more general approach, 
when there is a need to analyse phases which cannot be described by linear 
ODEs and how their limit distribution and the behavior of the trajectories 
will change. 
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Appendix A. Simplifying the system of ODEs 
and derivation of the formula for the attractor 


Let us consider the system ([^ with a given value of i, i = 0,1. Equivalently 
we can rewrite it as x' = Mx -|- c, where 


■ -1 

0 

0 


■ 1 ■ 

a 

—a 

0 

and c = 

0 

0 

b 

-b 


0 
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We notice that M has three distinct eigenvalues: —1, —a, —b and we can 
choose the eigenvectors vi, V 2 , V 3 in such a way that the vector 1 = [ 1 , 1 , 1 ] 
has also the coordinates [1,1,1] in the basis of the eigenvectors. Precisely, 
the eigenvectors vi, V 2 , V 3 are given by the formulas: 


Vi = 

V2 = 

V3 = 


1 , 


0 , 


ab 


a — 1 ’ (a — 1)(6 — 1 ) 
-1 b 


a — 1 ’ (a — l)(a — 6 ) 
ab 

0 , 0 , 1 - 


(a — 1)(6 — 1 ) (a 


b 

l)(a — b). 


( 22 ) 


We transform (0 by rewriting it in the new basis and we obtain new formulas 
for the system!^: 

(^^-ax2 + ai (23) 

for i = 0,1. Let x = [xi,X 2 ,X 3 \'^ be a column vector and let vrj(x) denote 
the solution of (|^ at time t with the initial condition x. We get 

7 r°(x) = (e“*Xi,e"“*X 2 ,e““x 3 ) (24) 


and 

TTt (x) = 1 + Tit (x) - 7 rt°(l), (25) 

where 1 denotes now a column vector [1,1,1]^. By alternate compositing of 
these functions we have the formulas 


K (x) = 1 + K+t2 (x) - 1 (26) 

as well 

(^) = 1 + ^l+t2 (^) - K+t2 1 ( 27 ) 

for any ^ 1,^2 > 0. Now, substituting := a, := (5 we get 


^iVo(x) = {a + I3{xi - 1),q;“ + /l“(a:2 - 1), a’’ + (3''{x3 - 1)), (28) 

and 

(x) = (1 - CK + f3xi, 1 - «“ + 13°'X 2 , 1 - CK^ + /S’^Xs) (29) 
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where 1 > a > /3 > 0. Taking as the initial points x = (0, 0, 0) in the formula 
(28) and x = (1,1,1) in the formula (29), we get parametric equations for 
the surfaces Aq and Ai which are the boundaries of A (see Sec. |^: 


Ao = {(a — /3, — /3“, — /3^) : 1 > a > /9 > 0}, 

Ai = {(1 — a + /9,1 — + /3“, 1 — + /3^) : 1 > a > /3 > 0}. 


Similarly: 


(x) = (a - /3 + 72^1, - /?“ + 7 x 2 , 

W = (1 - a + /3 + 7(xi - 1), 1 - + /3“ + 7“(a:2 - 1), 1 - + 7'’(a;3 - 1)), 

where l>a>/d> 7>0. 

Now, let Id = {{x,y, z): 0 < z < y < X < 1} and we define the function 
f; Id —)■ by the formula 

f (a:, y,z) = {x — y + z,x°' — y°' + 2:“, x^ — y’’ + z’’). (30) 

It is easy to check that f is a local diffeomorphism. In particular, 

f(i/) = {x = 7r° TT^^TT^ 1 : ti > 0, ta > 0, > 0} (31) 

is an open set and f(Id) is the interior of A. Hence 

A = {(x — y + z, a:“ — |/“ + 2:“, x^ — y^ + z^) ■. 1x y ^ z ^ Q} (32) 

and Aq i Ai are indeed the boundaries of A. Moreover Ai is the symmetrical 
image of Aq (and vice versa) with respect to a point (|, |, |) • From this 
property we have the equivalent dehnition of A, i.e.: 

A = {(l-(x-y+z), l-(a;“-2/“ + z“), 1 -(x’’-y’’+ z^) : 1 ^ x ^ y ^ z ^ 0}. 

(33) 

Appendix B. The proof of asymptotic stability 

Now we will prove the main result of this paper. We use the following theo¬ 
rem. 
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Theorem 2 Let 'K be a compact metric space and S be the Borel a—algebra. 
If a Markov semigroup {-P(t)}t^o satisfies two conditions: 

(a) for every density f we have JTmfdt > 0 a.e., 

(b) for every qq E X there exist k > 0, t > 0 and a measurable function 
rj ^ 0 such that f 7](p) m(dp) > 0 and 

^(i)f(p) ^ V{P) [ fiq)midq), 

J B(qo,K) 

for p G X, where B{qo, tf) is the open ball with center go o,nd radius k, 
then the semigroup {P{t)}t^Q is asymptotically stable. 

Theorem was formulated as Corollary 1 in and follows from two earlier 
results: 

Theorem 3 If {P{t)}t>o is a partially integral Markov semigroup and 
has a unique invariant density /* > 0, then the semigroup {P(f)}t>o is 
asymptotically stable. 

Theorem 4 If {P{t)}t>o is a Markov semigroup on a metric space, 

satisfies (a) and (b), and has no invariant density, then it is sweeping from 
compact sets. 

From Theorem it follows that if the space X is compact and the semigroup 
satisfies (a) and (b), then it has a unique and positive invariant density. Now 
Theorem follows immediately from Theorem 

Hence, the idea of the proof is as follows. First, we check that all tra¬ 
jectories enter the set A and this set is invariant. This allows us to reduce 
the proof of asymptotic stability only to the semigroup {P{t)}t>o restricted 
to the set A. Then we check that conditions (a) and (b) are satished on A. 
Since A is compact the semigroup {P{t)}t>o is asymptotically stable. 

Firstly, we introduce some necessary dehnitions. 

Definition 6 Let V (M) be the set of real smooth vector fields on the mani¬ 
fold M on and let C°°{M) denote the set of a real-valued smooth functions 
on V{M). A Lie bracket of two vector fields a,b E V{M) is a vector field 
given by the formula: 
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Definition 7 Let a PDMP be defined by the systems of differential equations 
x' = gfix), i E I = {0,1,/c}, k E N. We say that the Hormander’s 
condition holds at a point x if vectors 

92^^) ■ ■ ■ 1 \9i: 9j\i.^^l<i,j<k: \9ii \9 j i 9i\\^^^l<i,j,l<ki • • • 

span the space 

Definition 8 Let n G N, t > 0, r = (ri, T 2 , ■ ■ ■ ,t — Tn-i — ■■■ — Ti) and 
i = {ii,... ,in) such that for all k E {I, ..., n — 1} we have Tk > 0, 4 7^ ik+i 
and ik E {0,1}. A function 

f’xuAr) := o ttAz] o ■ ■ ■ o 

is called a cumulative flow along the trajectiories of the flows tt^L ..., tt*" 
with starting point x. 


Definition 9 We say that a point x E X communicates with y E X if there 
exist n G N, t > 0, r = (ti, r 2 ,..., t — r„_i — ■ ■ ■ — ti) and i = (zi ,... ,in) 
such that AxyAA — V- 

If every two points from the interior of X communicate, we call this property 
communication between states of the process. If for go E X there exists p E X 
such that go communicates with p and the Hormander’s condition holds at 
the point p, then go satisfies condition (b). This fact is a simple consequence 
of Theorem 4], 

Let us denote by aflx) a vector field representing the system (|^ with a 
fixed value of z G {0,1} at a point x G [0,1]^. After short calculation of the 
following expressions: 

ai-ao = (1, 0, 0), [ao, Oi] = (1, -a, 0), [oo, [oo, Oi]] = (1, -(a^ + a), ab); 

we obtain three linear independent vectors in space. Hence, these vectors 
span and condition (b) of Theorem]^ holds. However, it gets more difficult 
to check condition (a), because it does not hold on the whole space [0,1]^ x 
{0,1}. We will prove that (a) holds on A. Moreover, A is a stochastic 
attractor, i.e. a measurable subset of [0,1]^ such that for every density / G 
L^{A) we have 


lim 

t^OO 


F(f)/(x, z) dx di = lim P(C(^) G A) = 1. 


(34) 



First, we show that A is an invariant set for the process. It follows from 
the fact that if we take any x G A, then we stay in A under the action of 
both semi-flows given by (|^ with the initial condition x. In other words, we 
check that if we take any x G A and f > 0, then both 7r°(x) and (x) stay in 
A. Since the process switches between these two flows, its trajectories cannot 
leave A. 

Let us dehne the set: 


D = {{x—y+z—w, x^—y’^+z’^—w^): l^x^y^z^w^O}. 

We prove that D is the same set as A. For s G {1, a, 6} and any 

do = (a:o -yo + zo- Wo, x^ - y^ + z^ - x^-y^ + z^- Wg) G D 

we have 

- 2/0 + ^0 - ^ J " Vo + 4') 

\ -Zq -Zq -Zq / 

= (xo^o)* - (a^o2/o)" + i^oZoY, 

where the second equality follows from the equivalence of two dehnitions of 
the set A (see 33). Therefore, A = D. Since it is easy to check that for any 
f > 0 and i G {0,1} we have 7il{A) C A*, A is an invariant set for our process. 

there exist > 0, 2° > 0 such that we 


Using the formula 
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have 7r°o7r^\)7r°ol = (|, Y- From continuous dependence of the solutions on 
the initial condition, we can hnd 5 > 0 and e > 0 such that for every y G 
with ||y — 1|| <6 and t = (A, t 2 , A) with |A — < £ for i = 1, 2, 3 we have 

^ Int A. Moreover, there exists A > 0 such that ~ 1|| < 

6 for each point x G [0,1]^ and t > to- Thus 7rj(^(x) G IntA if 

X G [0,1]^ and |A — 2°] < £ for i = 1,2,3 and t > to- The probability 
that the sequence of the consecutive hve jump moments Tn,, Tn +4 has the 
properties: T„+i - > A and T^+i+i - T^+i e (t° - for * = 1,2,3; 

is bounded from below (see [^) by some positive number r/ > 0. Since our 
PDMP enters the cube [0,1]^ (see the discussion under the formula ([^) we 
have also that it enters the interior of the attractor with probability one. 


which completes the proof of (34). 


To prove the remaining condition (a), we use the fact that it is equivalent 
to communication between states for x, y G intA and hxed i G {0,1}. In 
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other words, we show that the cumulative flow consisting the flows of vr® 
and TT^ with the initial condition x G int A generates whole A. The question 
to face is: does the total control of the system (0 between two arbitrary 
points in the interior of A exist? The problem of control for linear dynamical 
systems has been extensively studied in the past ( 0. Ezl). but this special 
case appears to be relatively far from classical results of the controllability 
theory and seems to not undergo any of those procedures. However, the 
proof of the communication between states property in this special case is 
surprisingly simple. Due to symmetry of A, we consider only these cumulative 
flows which begin from 7r°. The case when we start from tt^ is analogous. 
Fix a point x G H. After compositing four transformations we obtain 


V’x,t,0 := TT^TT^J 


= 1 - TT.^l + + ^u+t,+ts+tA^) 


Fix e: > 0 and take any points x,y,z such that s<z<y<x<l. 
Then we can And ti > 0, for i = 1,2, 3,4, such that x = y = 
z = ^ Thus 


'0x,t,o = {I - X + y - z + exiA - x"" + y"" - z"" + e°'X 2 , l-x^ + y'' - z^ + 

Let 

A^ = {{1 — X + y — z,l — x°' + y^ — z'^ ,1 — x’^ + y^ — z’^): e<z<y<x<l} 

and Ve = {exi,e°‘X 2 ,s^X 3 ). Then, starting from the point x and using a 
composition of four transformations we communicate with each point from 
the set 

A, + v, := {y+ v,: y G A,}. 

Since Int Ac IJ A^, we conclude that we can join x and any interior point of 

£>0 

A by '0x,t,o (note that this property may not hold for the boundary points of 
A). As a result, condition (a) from Theorem]^ is satisfied and the semigroup 
{P(t)}t^o is asymptotically stable. 

Remark 3. The proof in the general case (i.e. including the subcases 
a = 6, a = 1 and 6 = 1) is similar to the presented above, but technically more 
difficult. We do not change the variables in the system ([^ (see Appendix 
A), but we define the attractor for the process as the closure of the set: 

Al = {x = TT^jl - TT^^l + TT^gl : 0 < Si < S2 < S3} ( 35 ) 
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or equivalently as the closure of the set: 


^2 = {x = 1 - 7r°^l + TT^^l - 7r°3l : 0 < Si < Sa < S 3 }. (36) 

We have clAi = CIA 2 , because both sets Ai and A 2 have the same bound¬ 
aries: 

cl {x = 1 - 7 r° 3 1 : 0 < Si < Sa} U cl {x = 1 - 7r°^l-h tt^^I : 0 < Si < Sa}. 

Having these formulas for the attractor one can prove that indeed this set 
is an attractor and that any two points in the interior of the attractor can 
communicate. 
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